#pragma once
#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <stdio.h>
#include <limits>
#include <array>
#include <cstddef>
#include <omp.h>
// #include <mpi.h>

template <typename V>
struct GridPoint
{
    /* data */
    V _x, _y, _z;
    int _interfaceIdx;

    GridPoint() : _x(0), _y(0), _z(0), _interfaceIdx(0) {}

    GridPoint(V xValue, V yValue, V zValue, int interfaceIdx = 0)
            : _x(xValue), _y(yValue), _z(zValue), _interfaceIdx(interfaceIdx) {}
    
    // 输出运算符
    friend std::ostream& operator<<(std::ostream& os, const GridPoint& point) {
        os << "(" << point._x << ", " << point._y << ", " << point._z << ", interfaceIdx: " << point._interfaceIdx << ")";
        return os;
    }
};

template <typename V>
V distance2(const GridPoint<V>& p1, const GridPoint<V>& p2)
{
    V dist2 = (p1._x - p2._x) * (p1._x - p2._x) + (p1._y - p2._y) * (p1._y - p2._y) + (p1._z - p2._z) * (p1._z - p2._z);
    return dist2;
}
template <typename V, typename T, std::size_t N>
struct GridCell
{
    std::array<V, N> _points;      // 网格顶点
    T _youngMod;                   // 杨氏模量
    T _poissonRatio;               // 泊松比
    T _rockDensity;                // 岩石密度
    T _porosity;                   // 孔隙率
    int _cellZone;                 // 单元区域索引

    // 默认构造函数
    GridCell()
        : _points(), _youngMod(0.0), _poissonRatio(0.0), _rockDensity(0.0), _porosity(0.0), _cellZone(0) {}

    // 带参数的构造函数
    GridCell(const std::array<V, N>& points, T youngMod = 0.0, T poissonRatio = 0.0, T rockDensity = 0.0, T porosity = 0.0, int cellZone = 1)
        : _points(points), _youngMod(youngMod), _poissonRatio(poissonRatio), _rockDensity(rockDensity), _porosity(porosity), _cellZone(cellZone) {}

    GridCell(const GridCell& other) = default;

    GridCell& operator=(const GridCell& other) = default;

    ~GridCell() = default;

    // 输出流运算符
    friend std::ostream& operator<<(std::ostream& os, const GridCell& cell) {
        os << "GridCell(points: [";
        for (size_t i = 0; i < N; ++i) {
            os << cell._points[i];
            if (i < N - 1) os << ", ";
        }
        os << "], youngMod: " << cell._youngMod
           << ", poissonRatio: " << cell._poissonRatio
           << ", rockDensity: " << cell._rockDensity
           << ", porosity: " << cell._porosity
           << ", cellZone: " << cell._cellZone << ")";
        return os;
    }
};

template <typename V>
struct Fault
{   
    std::string _name; // 断层名称
    std::vector<std::vector<GridPoint<V>>> _lines; // 断层柱线

    Fault() = default;

    void add_line(const std::vector<GridPoint<V>>& faultLine) {
        _lines.push_back(faultLine);
    }

    // 获取当前柱线数
    size_t line_count() const {
        return _lines.size();
    }
};

/**
 * @class CornerGrid
 * @brief A class for handling corner-point grids used in geological modeling and simulation.
 *
 * This class is designed to parse, construct, and manipulate corner-point grid data. 
 * It supports grid structure assembly, property reading, fault integration, 
 * and exporting grid data to VTK format for visualization.
 *
 * Template Parameters:
 * - `CellType` : The type of cell structure in the grid (e.g., `GridCell`).
 * - `DataType` : The type of numerical values used for grid points and properties (e.g., `float` or `double`).
 *
 * Features:
 * - Reads grid structure and properties from industry-standard file formats.
 * - Supports the integration of geological faults into the grid structure.
 * - Exports the grid and its properties to VTK format for 3D visualization.
 * - Provides utility functions for querying grid dimensions and faults.
 *
 * Example Usage:
 * ```
 * CornerGrid<GridCell<GridPoint<double>, double, 8>, double> grid("GridName");
 * grid.LoadGrid("grid_file.out");
 * grid.LoadGridProperties("property_file.txt");
 * grid.LoadFault("fault_file.txt");
 * ```
 */
template <bool IsFineGrid, typename DataType, typename CellType>
class CornerGrid
{   
private:
    std::string _gridName;
    int _Nx, _Ny, _Nz;
    DataType _xOffset, _yOffset;
    DataType _xmax, _xmin, _ymax, _ymin, _zmax, _zmin;
    std::vector<CellType> _cellArray;
    std::vector<Fault<DataType>> _faultArray;
    // std::vector<grid_point<T>> pointarray;

    // Function to read the grid size from a file (SPECGRID keyword section)
    // NX, NY, NZ represent the dimensions of the grid
    void _read_gridinfo(const std::string filename, int& NX, int& NY, int& NZ, DataType& X_offset, DataType& Y_offset)
    {
        std::ifstream file(filename);
        if (!file.is_open()) // Check if the file can be opened
        {
            std::cerr << "Error: could not open file " << filename << std::endl;
            return;
        }
        std::string line;
        while (std::getline(file, line)) // Iterate through each line of the file
        {   
            if (line.find("MAPAXES") != std::string::npos)
            {
                std::getline(file, line);
                std::getline(file, line);

                std::istringstream iss(line);
                iss >> X_offset >> Y_offset;
            }
            if (line.find("SPECGRID") != std::string::npos) // Find the SPECGRID section
            {
                std::getline(file, line); // Read the next line containing grid dimensions
                std::istringstream iss(line);
                iss >> NX >> NY >> NZ; // Extract grid dimensions
                break;
            }
        }
    }

    void _read_gridpoints(const std::string filename, const int NX, const int NY, const int NZ, 
                            std::vector<GridPoint<DataType>>& coords, std::vector<DataType>& zcorns)
    {
        std::ifstream file(filename);
        if (!file.is_open()) // Check if the file can be opened
        {
            std::cerr << "Error: could not open file " << filename << std::endl;
            return;
        }
        std::string line;

        // Read COORD data (X, Y coordinates of grid corners)
        while (std::getline(file, line))
        {
            if (line.find("COORD") != std::string::npos)
            {
                for (int i = 0; i < (NX + 1) * (NY + 1); i++)
                {   
                    DataType x1, y1, z1, x2, y2, z2;

                    std::getline(file, line);
                    std::istringstream iss(line);
                    iss >> x1 >> y1 >> z1 >> x2 >> y2 >> z2;
                    coords[i * 2] = GridPoint<DataType>(x1, y1, z1);
                    coords[i * 2 + 1] = GridPoint<DataType>(x2, y2, z2);
                }
                break;
            }
        }
        
        // Read ZCORN data (Z coordinates of grid corners)
        DataType z;
        int zcorn_count = 0;
        while (std::getline(file, line))
        {
            if (line.find("ZCORN") != std::string::npos)
            {
                std::getline(file, line, '/'); // Extract data up to delimiter '/'
                std::istringstream iss(line);
                while (iss >> z)
                {
                    zcorns[zcorn_count++] = z; // Store Z coordinate
                }
                break;
            }   
        }

        int expected_zcorn_size = NX * NY * NZ * 8;
        if (zcorn_count != expected_zcorn_size)
        {
            std::cerr << "Error: ZCORN data size mismatch. Expected " << expected_zcorn_size 
                    << " but got " << zcorn_count << std::endl;
        }
    }

    void _read_propdata(const std::string filename, std::vector<CellType>& cells)
    {
        std::ifstream file(filename);
        if (!file.is_open()) // Check if the file can be opened
        {
            std::cerr << "Error: could not open file " << filename << std::endl;
            return;
        }

        // Declare a string to store lines read from the file.
        std::string line;

        // Search for the section containing YOUNGMOD
        while (std::getline(file, line))
        {
            if (line.find("YOUNGMOD") != std::string::npos)
            {
                // Read the line of data after the "YOUNGMOD" keyword until the delimiter '/'.
                std::getline(file, line, '/');
                std::istringstream iss(line);  // Create a string stream for parsing the data.
                std::string tmp;              // Temporary string to hold each token.
                int cnt = 0;                  // Counter to track the position in the `youngmods` array.

                // Parse the data in the format `<num>*<value>` (e.g., `3*5.2`).
                while (iss >> tmp)
                {
                    // Split the token into `<num>` and `<value>` parts.
                    auto pos = tmp.find('*');           // Find the position of the '*' character.
                    int num = std::stoi(tmp.substr(0, pos)); // Extract `<num>` (number of repetitions).
                    DataType val = std::stof(tmp.substr(pos + 1));  // Extract `<value>` (actual data value).

                    // Fill the `youngmods` array with `<value>` repeated `<num>` times.
                    for (int i = 0; i < num; i++)
                    {
                        cells[cnt++]._youngMod = val;  // Assign the value to the array and increment the counter.
                    }
                }
            }

            if (line.find("POISSONR") != std::string::npos)
            {
                // Read the line of data after the "POISSONR" keyword until the delimiter '/'.
                std::getline(file, line, '/');
                std::istringstream iss(line);  // Create a string stream for parsing the data.
                std::string tmp;              // Temporary string to hold each token.
                int cnt = 0;                  // Counter to track the position in the `poissonrs` array.

                while (iss >> tmp)
                {
                    // Split the token into `<num>` and `<value>` parts.
                    auto pos = tmp.find('*');           // Find the position of the '*' character.
                    int num = std::stoi(tmp.substr(0, pos)); // Extract `<num>` (number of repetitions).
                    DataType val = std::stof(tmp.substr(pos + 1));  // Extract `<value>` (actual data value).

                    // Fill the `poissonrs` array with `<value>` repeated `<num>` times.
                    for (int i = 0; i < num; i++)
                    {
                        cells[cnt++]._poissonRatio = val;  // Assign the value to the array and increment the counter.
                    }
                }
            }

            if (line.find("ROCKDEN") != std::string::npos)
            {
                std::getline(file, line, '/');
                std::istringstream iss(line);
                std::string tmp;
                int cnt = 0;

                while (iss >> tmp)
                {
                    auto pos = tmp.find('*');
                    int num = std::stoi(tmp.substr(0, pos));
                    DataType val = std::stof(tmp.substr(pos + 1));

                    for (int i = 0; i < num; i++)
                    {
                        cells[cnt++]._rockDensity = val;
                    }
                }
            }

            if (line.find("PORO") != std::string::npos)
            {
                std::getline(file, line, '/');
                std::istringstream iss(line);
                std::string tmp;
                int cnt = 0;

                while (iss >> tmp)
                {
                    auto pos = tmp.find('*');
                    int num = std::stoi(tmp.substr(0, pos));
                    DataType val = std::stof(tmp.substr(pos + 1));

                    for (int i = 0; i < num; i++)
                    {
                        cells[cnt++]._porosity = val;
                    }
                }
            }
        }
    }

    void _read_fault(const std::string filename, std::vector<Fault<DataType>>& faults, DataType X_offset = 0, DataType Y_offset = 0)
    {
        std::ifstream file(filename);
        if (!file.is_open()) // Check if the file can be opened
        {
            std::cerr << "Error: could not open file " << filename << std::endl;
            return;
        }

        std::string line;       // Declare a string to store lines read from the file.
        Fault<DataType>* current_fault = nullptr;
        std::vector<GridPoint<DataType>> current_line;
        int current_line_number = -1;

        while (std::getline(file, line))
        {   
            if (line.empty() || line.find_first_not_of(" \t") == std::string::npos) {
                continue;
            }
            std::istringstream iss(line);

            if (line.find("FAULT") == 0) {
                std::string keyword, fault_name;
                iss >> keyword >> fault_name; // 读取 "FAULT" 和断层名称

                // 将上一断层的最后一条柱线存入
                if (current_fault && !current_line.empty()) {
                    current_fault->add_line(current_line);
                    current_line.clear();
                }

                faults.emplace_back();
                current_fault = &faults.back();
                current_fault->_name = fault_name;

                current_line_number = -1;
            } else if (current_fault) {
                // 解析点数据行
                DataType x, y, z;
                int line_number;
                iss >> x >> y >> line_number >> z;

                if (line_number != current_line_number) {
                    if (!current_line.empty()) {
                        current_fault->add_line(current_line); 
                        current_line.clear();
                    }
                    current_line_number = line_number;
                }

                // 将点加入当前断层的对应柱线
                current_line.push_back(GridPoint<DataType>(x - X_offset, Y_offset - y, -z));
            }
        }

        if (current_fault && !current_line.empty()) {
            current_fault->add_line(current_line);
        }
    }

    // void _assemble_points(const std::vector<grid_point<T>>& coords, const std::vector<T>& zcorns, int NX, int NY, int NZ, std::vector<grid_point<T>>& points)
    // {   
    //     for (int k = 0; k < NZ; k++)
    //     {
    //         for (int j = 0; j < NY; j++)
    //         {
    //             for (int i = 0; i < NX; i++)
    //             {
    //                 int cellIndex = k * NY * NX + j * NX + i; // Compute cell index
    //                 for (int index = 0; index < 8; index++) // Process 8 corners of the cell
    //                 {
    //                     int x = i + (index & 1); // Determine X offset
    //                     int y = j + ((index & 2) >> 1); // Determine Y offset
    //                     grid_point<T> p1 = coords[(y * (NX + 1) + x) * 2];
    //                     grid_point<T> p2 = coords[(y * (NX + 1) + x) * 2 + 1];

    //                     T pz = zcorns[k * NY * NX * 8 + ((index & 4) >> 2) * NX * NY * 4 + j * NX * 4 + ((index & 2) >> 1) * NX * 2 + i * 2 + (index & 1)];
    //                     T px = p1.x + (pz - p1.z) * (p2.x - p1.x) / (p2.z - p1.z);
    //                     T py = p1.y + (pz - p1.z) * (p2.y - p1.y) / (p2.z - p1.z);
    //                     points[cellIndex * 8 + index] = grid_point<T>(px, py, pz);
    //                 }
    //             }
    //         }
    //     }
    // }

    void _assemble_cells(const std::vector<GridPoint<DataType>>& coords, const std::vector<DataType>& zcorns, int NX, int NY, int NZ, std::vector<CellType>& cells)
    {
        #pragma omp parallel for collapse(1)
        for (int k = 0; k < NZ; k++) {
            int interface_idx = 11 + NZ - k, cellzone_idx;
            if (k == 0) { cellzone_idx = 5; }
            else if (k == NZ - 1) { cellzone_idx = 3; }
            else { cellzone_idx = 1; }

            for (int j = 0; j < NY; j++) {
                for (int i = 0; i < NX; i++) {
                    int cellIndex = k * NY * NX + j * NX + i; // Compute cell index
                    for (int index = 0; index < 8; index++) // Process 8 corners of the cell
                    {
                        int x = i + (index & 1); // Determine X offset
                        int y = j + ((index & 2) >> 1); // Determine Y offset
                        GridPoint<DataType> p1 = coords[(y * (NX + 1) + x) * 2];
                        GridPoint<DataType> p2 = coords[(y * (NX + 1) + x) * 2 + 1];

                        DataType pz = zcorns[k * NY * NX * 8 + ((index & 4) >> 2) * NX * NY * 4 + j * NX * 4 + ((index & 2) >> 1) * NX * 2 + i * 2 + (index & 1)];
                        DataType px = p1._x + (pz - p1._z) * (p2._x - p1._x) / (p2._z - p1._z);
                        DataType py = p1._y + (pz - p1._z) * (p2._y - p1._y) / (p2._z - p1._z);
                        if (px < _xmin) _xmin = px;
                        if (px > _xmax) _xmax = px;
                        if (py < _ymin) _ymin = py;
                        if (py > _ymax) _ymax = py;
                        if (pz < _zmin) _zmin = pz;
                        if (pz > _zmax) _zmax = pz;
                        cells[cellIndex]._points[index] = GridPoint<DataType>(px, py, pz, interface_idx - ((index & 4) >> 2));
                        // points[cellIndex * 8 + index] = GridPoint<T>(px, py, pz);
                    }
                    cells[cellIndex]._cellZone = cellzone_idx;
                }
            }
        }
    }
    
    void _write_vtk(const std::string filename, const std::vector<CellType>& cells)
    {
        FILE* file = fopen(filename.c_str(), "w");
        if (file == NULL)
        {
            std::cerr << "Error: could not open file " << filename << std::endl;
            return;
        }

        fprintf(file, "# vtk DataFile Version 3.0\n");
        fprintf(file, "%s in vtk\n", _gridName.c_str());
        fprintf(file, "ASCII\n");
        fprintf(file, "DATASET UNSTRUCTURED_GRID\n\n");

        int totalPoints = static_cast<int>(cells.size()) * 8;

        fprintf(file, "POINTS %d float\n", totalPoints);

        // 按顺序写入所有点的坐标
        for (size_t i = 0; i < cells.size(); ++i) {
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[0]._x), static_cast<DataType>(cells[i]._points[0]._y), static_cast<DataType>(cells[i]._points[0]._z));
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[1]._x), static_cast<DataType>(cells[i]._points[1]._y), static_cast<DataType>(cells[i]._points[1]._z));
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[2]._x), static_cast<DataType>(cells[i]._points[2]._y), static_cast<DataType>(cells[i]._points[2]._z));
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[3]._x), static_cast<DataType>(cells[i]._points[3]._y), static_cast<DataType>(cells[i]._points[3]._z));
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[4]._x), static_cast<DataType>(cells[i]._points[4]._y), static_cast<DataType>(cells[i]._points[4]._z));
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[5]._x), static_cast<DataType>(cells[i]._points[5]._y), static_cast<DataType>(cells[i]._points[5]._z));
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[6]._x), static_cast<DataType>(cells[i]._points[6]._y), static_cast<DataType>(cells[i]._points[6]._z));
            fprintf(file, "%f %f %f\n", static_cast<DataType>(cells[i]._points[7]._x), static_cast<DataType>(cells[i]._points[7]._y), static_cast<DataType>(cells[i]._points[7]._z));
        }

        int numCells = static_cast<int>(cells.size());
        fprintf(file, "\nCELLS %d %d\n", numCells, numCells * 9);  // 每个单元包含1个单元点数量信息和8个点ID，总共9个整数

        // 写入每个单元的点索引
        int pointIdOffset = 0;
        for (int i = 0; i < numCells; ++i) {
            fprintf(file, "8");  // 8个点的单元（六面体）
            fprintf(file, " %d", pointIdOffset + 0);
            fprintf(file, " %d", pointIdOffset + 1);
            fprintf(file, " %d", pointIdOffset + 2);
            fprintf(file, " %d", pointIdOffset + 3);
            fprintf(file, " %d", pointIdOffset + 4);
            fprintf(file, " %d", pointIdOffset + 5);
            fprintf(file, " %d", pointIdOffset + 6);
            fprintf(file, " %d", pointIdOffset + 7);
            fprintf(file, "\n");
            pointIdOffset += 8;
        }

        fprintf(file, "\nCELL_TYPES %d\n", numCells);
        for (int i = 0; i < numCells; ++i)
        {
            fprintf(file, "11\n");
        }

        fprintf(file, "\nPOINT_DATA %d\n", totalPoints);
        fprintf(file, "SCALARS IsInterface int 1\n");
        fprintf(file, "LOOKUP_TABLE default\n");
        // TODO: Write interface data
        for (size_t i = 0; i < cells.size(); ++i) {
            fprintf(file, "%d\n", cells[i]._points[0]._interfaceIdx);
            fprintf(file, "%d\n", cells[i]._points[1]._interfaceIdx);
            fprintf(file, "%d\n", cells[i]._points[2]._interfaceIdx);
            fprintf(file, "%d\n", cells[i]._points[3]._interfaceIdx);
            fprintf(file, "%d\n", cells[i]._points[4]._interfaceIdx);
            fprintf(file, "%d\n", cells[i]._points[5]._interfaceIdx);
            fprintf(file, "%d\n", cells[i]._points[6]._interfaceIdx);
            fprintf(file, "%d\n", cells[i]._points[7]._interfaceIdx);
        }

        fprintf(file, "\nCELL_DATA %d\n", numCells);

        fprintf(file, "SCALARS CellZone int 1\n");
        fprintf(file, "LOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i)
        {
            fprintf(file, "%d\n", cells[i]._cellZone);
        }

        // 写入单元数据（例如YOUNGMOD, POISSONR, ROCKDEN, PORO）
        // if constexpr (std::is_same_v<CellType, GridCell<GridPoint<DataType>, DataType, 8>>)
        if constexpr (IsFineGrid)
        {
            fprintf(file, "SCALARS YoungModulos float\n");
            fprintf(file, "LOOKUP_TABLE default\n");
            for (int i = 0; i < numCells; ++i)
            {
                fprintf(file, "%f\n", static_cast<DataType>(cells[i]._youngMod));
            }

            // 写入POISSONR属性
            fprintf(file, "\nSCALARS PoissonRatio float\n");
            fprintf(file, "LOOKUP_TABLE default\n");
            for (int i = 0; i < numCells; ++i)
            {
                fprintf(file, "%f\n", static_cast<DataType>(cells[i]._poissonRatio));
            }

            // 写入ROCKDEN属性
            fprintf(file, "\nSCALARS RockDensity float\n");
            fprintf(file, "LOOKUP_TABLE default\n");
            for (int i = 0; i < numCells; ++i)
            {
                fprintf(file, "%f\n", static_cast<DataType>(cells[i]._rockDensity));
            }

            // 写入PORO属性
            fprintf(file, "\nSCALARS Porosity float\n");
            fprintf(file, "LOOKUP_TABLE default\n");
            for (int i = 0; i < numCells; ++i)
            {
                fprintf(file, "%f\n", static_cast<DataType>(cells[i]._porosity));
            }
        }

        fclose(file);
    }

    void _write_fault_vtk(const std::string filename, const std::vector<Fault<DataType>>& faults)
    {
        FILE* file = fopen(filename.c_str(), "w");
        if (file == NULL)
        {
            std::cerr << "Error: could not open file " << filename << std::endl;
            return;
        }

        fprintf(file, "# vtk DataFile Version 3.0\n");
        fprintf(file, "%s_Fault in vtk\n", _gridName.c_str());
        fprintf(file, "ASCII\n");
        fprintf(file, "DATASET UNSTRUCTURED_GRID\n\n");

        // 收集所有点并分配索引
        std::vector<GridPoint<DataType>> all_points;
        std::vector<std::vector<std::vector<size_t>>> fault_line_point_indices;

        for (const auto& fault : faults) {
            std::vector<std::vector<size_t>> fault_point_indices;

            for (const auto& line : fault._lines) {
                std::vector<size_t> line_point_indices;
                for (const auto& point : line) {
                    all_points.push_back(point);
                    line_point_indices.push_back(all_points.size() - 1);
                }
                fault_point_indices.push_back(line_point_indices);
            }

            fault_line_point_indices.push_back(fault_point_indices);
        }

        // 写入 POINTS 部分
        fprintf(file, "POINTS %zu float\n", all_points.size());
        for (const auto& point : all_points) {
            fprintf(file, "%f %f %f\n", point._x, point._y, point._z);
        }

        // 计算单元数量和总大小
        size_t num_line_cells = 0;
        size_t total_line_size = 0;
        size_t num_surface_cells = 0;
        size_t total_surface_size = 0;

        // 计算线单元和面单元的数量和大小
        for (const auto& fault_point_indices : fault_line_point_indices) {
            // 线单元
            for (const auto& line_point_indices : fault_point_indices) {
                num_line_cells += 1;
                total_line_size += (line_point_indices.size() + 1);
            }
            // 面单元
            size_t num_lines_in_fault = fault_point_indices.size();
            if (num_lines_in_fault >= 2) {
                for (size_t i = 0; i < num_lines_in_fault -1; ++i) {
                    const auto& line1 = fault_point_indices[i];
                    const auto& line2 = fault_point_indices[i +1];

                    size_t n_quads = std::min(line1.size(), line2.size()) - 1;
                    for (size_t j = 0; j < n_quads; ++j) {
                        num_surface_cells += 1;
                        total_surface_size += (4 + 1);
                    }
                }
            }
        }

        size_t total_cells = num_line_cells + num_surface_cells;
        size_t total_size = total_line_size + total_surface_size;

        fprintf(file, "CELLS %zu %zu\n", total_cells, total_size);

        // 写入线单元
        for (const auto& fault_point_indices : fault_line_point_indices) {
            for (const auto& line_point_indices : fault_point_indices) {
                fprintf(file, "%zu", line_point_indices.size());
                for (const auto& idx : line_point_indices) {
                    fprintf(file, " %zu", idx);
                }
                fprintf(file, "\n");
            }
        }

        // 写入面单元（四边形）
        for (const auto& fault_point_indices : fault_line_point_indices) {
            size_t num_lines_in_fault = fault_point_indices.size();
            if (num_lines_in_fault >= 2) {
                for (size_t i = 0; i < num_lines_in_fault -1; ++i) {
                    const auto& line1 = fault_point_indices[i];
                    const auto& line2 = fault_point_indices[i +1];

                    size_t n_quads = std::min(line1.size(), line2.size()) -1;
                    for (size_t j = 0; j < n_quads; ++j) {
                        size_t p0 = line1[j];
                        size_t p1 = line1[j+1];
                        size_t p2 = line2[j+1];
                        size_t p3 = line2[j];

                        fprintf(file, "4 %zu %zu %zu %zu\n", p0, p1, p2, p3);
                    }
                }
            }
        }

        // 写入 CELL_TYPES 部分
        fprintf(file, "CELL_TYPES %zu\n", total_cells);

        // 线单元的类型代码为 4（VTK_POLY_LINE）
        for (size_t i = 0; i < num_line_cells; ++i) {
            fprintf(file, "4\n");
        }

        // 面单元的类型代码为 9（VTK_QUAD）
        for (size_t i = 0; i < num_surface_cells; ++i) {
            fprintf(file, "9\n");
        }

        fclose(file);
    }

public:
    CornerGrid(const std::string GridName) {
        _gridName = GridName;
    }

    void LoadGrid(const std::string gridFile, const std::string& propFile = "") {
        _xmax = -1e10, _xmin = 1e10;
        _ymax = -1e10, _ymin = 1e10;
        _zmax = -1e10, _zmin = 1e10;

        _read_gridinfo(gridFile, _Nx, _Ny, _Nz, _xOffset, _yOffset);

        std::vector<GridPoint<DataType>> coords((_Nx + 1) * (_Ny + 1) * 2);
        std::vector<DataType> zcorns(_Nx * _Ny * _Nz * 8);

        _read_gridpoints(gridFile, _Nx, _Ny, _Nz, coords, zcorns);

        // pointarray.resize(Nx * Ny * Nz * 8);
        _cellArray.resize(_Nx * _Ny * _Nz);

        // _assemble_points(coords, zcorns, Nx, Ny, Nz, pointarray);
        _assemble_cells(coords, zcorns, _Nx, _Ny, _Nz, _cellArray);

        if constexpr (IsFineGrid) {
            static_assert(std::is_same_v<CellType, GridCell<GridPoint<DataType>, DataType, 8>>, 
                          "LoadGrid is only available for GridCell types when loading properties!");
            if (!propFile.empty()) {
                _read_propdata(propFile, _cellArray);
            }           
        }
    }

    // void LoadGridProperties(const std::string prop_file) {
    //     static_assert(std::is_same_v<V, finegrid_cell<grid_point<T>, T, 8>>, 
    //                   "LoadGridProperties is only available for finegrid_cell types!");
    //     _read_propdata(prop_file, cellarray);
    // }

    void LoadFault(const std::string faultFile) {
        _read_fault(faultFile, _faultArray, _xOffset, _yOffset);
    }   

    void WriteToVtk(const std::string outputFile) {
        _write_vtk(outputFile, _cellArray);
    }

    void WriteFaultToVtk(const std::string outputFile) {
        _write_fault_vtk(outputFile, _faultArray);
    }    
    
    void add_cell(const CellType& cell) {
        _cellArray.push_back(cell);
    }

    void SetNx(int Nx) { _Nx = Nx; }
    void SetNy(int Ny) { _Ny = Ny; }
    void SetNz(int Nz) { _Nz = Nz; }
    void SetXOffset(DataType xOffset) { _xOffset = xOffset; }
    void SetYOffset(DataType yOffset) { _yOffset = yOffset; }

    int GetNx() const { return _Nx; }
    int GetNy() const { return _Ny; }
    int GetNz() const { return _Nz; }
    DataType GetXmax() const { return _xmax; }
    DataType GetXmin() const { return _xmin; }
    DataType GetYmax() const { return _ymax; }
    DataType GetYmin() const { return _ymin; }
    DataType GetZmax() const { return _zmax; }
    DataType GetZmin() const { return _zmin; }  
    DataType GetXOffset() const { return _xOffset; }
    DataType GetYOffset() const { return _yOffset; }
    std::string GetGridName() const { return _gridName; }
    std::vector<CellType> GetCellArray() const { return _cellArray; }
    std::vector<Fault<DataType>> GetFaultArray() const { return _faultArray; }
};

template <typename CellType, typename DataType, std::size_t N>
class FinElemGrid
{
private:
    std::string _gridName;
    std::vector<GridPoint<DataType>> _pointArray;
    std::vector<CellType> _cellArray;
    int _Num_OriginCell;
    // int _Overburden_offset, _Undersides_offset, _Underburden_offset, _Sidesburden_offset;

    void _write_vtk(const std::string filename, const std::vector<GridPoint<DataType>> points, const std::vector<CellType>& cells)
    {
        FILE* file = fopen(filename.c_str(), "w");
        if (file == NULL)
        {
            std::cerr << "Error: could not open file " << filename << std::endl;
            return;
        }

        fprintf(file, "# vtk DataFile Version 3.0\n");
        fprintf(file, "%s in vtk\n", _gridName.c_str());
        fprintf(file, "ASCII\n");
        fprintf(file, "DATASET UNSTRUCTURED_GRID\n\n");

        int totalPoints = static_cast<int>(points.size());
        fprintf(file, "POINTS %d float\n", totalPoints); 
        for (size_t i = 0; i < totalPoints; ++i) {
            const auto& pt = points[i];
            fprintf(file, "%f %f %f\n", static_cast<DataType>(pt._x), static_cast<DataType>(pt._y), static_cast<DataType>(pt._z));
        }

        int numCells = static_cast<int>(cells.size());
        fprintf(file, "\nCELLS %d %lu\n", numCells, numCells * (1 + N));  // 每个单元包含1个单元点数量信息和N个点ID，总共N个整数

        for (int i = 0; i < numCells; ++i) {
            fprintf(file, "%lu", N);
            for (size_t j = 0; j < N; ++j){
                fprintf(file, " %d", cells[i]._points[j]);
            }
            fprintf(file, "\n");
        }

        fprintf(file, "\nCELL_TYPES %d\n", numCells);
        for (int i = 0; i < numCells; ++i)
        {
            fprintf(file, "10\n");
        }

        fprintf(file, "\nPOINT_DATA %d\n", totalPoints);
        fprintf(file, "SCALARS IsInterface int 1\n");
        fprintf(file, "LOOKUP_TABLE default\n");

        for (size_t i = 0; i < totalPoints; ++i) {
            fprintf(file, "%d\n", points[i]._interfaceIdx);
        }

        fprintf(file, "\nCELL_DATA %d\n", numCells);
        fprintf(file, "SCALARS CellZone int 1\n");
        fprintf(file, "LOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i)
        {
            fprintf(file, "%d\n", cells[i]._cellZone);
        }

        fprintf(file, "SCALARS YoungModulos float\n");
        fprintf(file, "LOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i)
        {
            fprintf(file, "%f\n", static_cast<DataType>(cells[i]._youngMod));
        }

        // 写入POISSONR属性
        fprintf(file, "\nSCALARS PoissonRatio float\n");
        fprintf(file, "LOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i)
        {
            fprintf(file, "%f\n", static_cast<DataType>(cells[i]._poissonRatio));
        }

        // 写入ROCKDEN属性
        fprintf(file, "\nSCALARS RockDensity float\n");
        fprintf(file, "LOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i)
        {
            fprintf(file, "%f\n", static_cast<DataType>(cells[i]._rockDensity));
        }

        // 写入PORO属性
        fprintf(file, "\nSCALARS Porosity float\n");
        fprintf(file, "LOOKUP_TABLE default\n");
        for (int i = 0; i < numCells; ++i) {
            fprintf(file, "%f\n", static_cast<DataType>(cells[i]._porosity));
        }
        fclose(file);
    }

    // void _hexahedron26tetrahedron(const GridCell<GridPoint<DataType>, DataType, 8>& cell, size_t* tet_par)
    // {
    //     if (distance2(cell.points[1], cell.points[2]) >= distance2(cell.points[0], cell.points[3]))
    //     {
    //         tet_par[0] = 2;
    //         tet_par[1] = 3;
    //         tet_par[2] = 0;
    //         tet_par[3] = 4;

    //         tet_par[4] = 2;
    //         tet_par[5] = 4;
    //         tet_par[6] = 6;
    //         tet_par[7] = 3;

    //         tet_par[8] = 3;
    //         tet_par[9] = 4;
    //         tet_par[10] = 6;
    //         tet_par[11] = 7;

    //         tet_par[12] = 3;
    //         tet_par[13] = 5;
    //         tet_par[14] = 4;
    //         tet_par[15] = 7;

    //         tet_par[16] = 1;
    //         tet_par[17] = 3;
    //         tet_par[18] = 5;
    //         tet_par[19] = 4;

    //         tet_par[20] = 1;
    //         tet_par[21] = 0;
    //         tet_par[22] = 3;
    //         tet_par[23] = 4;
    //     }
    //     else
    //     {
    //         tet_par[0] = 2;
    //         tet_par[1] = 1;
    //         tet_par[2] = 0;
    //         tet_par[3] = 6;

    //         tet_par[4] = 0;
    //         tet_par[5] = 4;
    //         tet_par[6] = 6;
    //         tet_par[7] = 1;

    //         tet_par[8] = 1;
    //         tet_par[9] = 4;
    //         tet_par[10] = 6;
    //         tet_par[11] = 5;

    //         tet_par[12] = 1;
    //         tet_par[13] = 5;
    //         tet_par[14] = 6;
    //         tet_par[15] = 7;

    //         tet_par[16] = 1;
    //         tet_par[17] = 3;
    //         tet_par[18] = 7;
    //         tet_par[19] = 6;

    //         tet_par[20] = 1;
    //         tet_par[21] = 2;
    //         tet_par[22] = 3;
    //         tet_par[23] = 6;
    //     }
    // }

public:
    FinElemGrid(const std::string GridName) {
        _gridName = GridName;
    }

    void add_point(const GridPoint<DataType>& point) {
        _pointArray.push_back(point);
    }

    void add_cell(const CellType& cell) {
        _cellArray.push_back(cell);
    }

    void WriteToVtk(const std::string outputFile) {
        _write_vtk(outputFile, _pointArray, _cellArray);
    }

    void SetPointArray(const std::vector<GridPoint<DataType>>& pointArray) {
        _pointArray = pointArray;
    }
    void SetCellArray(const std::vector<CellType>& cellArray) {
        _cellArray = cellArray;
    }

    void SetNumOriginCell(int numOriginCell) { _Num_OriginCell = numOriginCell; }
    // void SetOverburdenOffset(int overburden_offset) { _Overburden_offset = overburden_offset; }
    // void SetUndersidesOffset(int undersides_offset) { _Undersides_offset = undersides_offset; }
    // void SetUnderburdenOffset(int underburden_offset) { _Underburden_offset = underburden_offset; }
    // void SetSidesburdenOffset(int sidesburden_offset) { _Sidesburden_offset = sidesburden_offset; }

    int GetNumOriginCell() const { return _Num_OriginCell; }
    // int GetOverburdenOffset() const { return _Overburden_offset; }
    // int GetUndersidesOffset() const { return _Undersides_offset; }
    // int GetUnderburdenOffset() const { return _Underburden_offset; }
    // int GetSidesburdenOffset() const { return _Sidesburden_offset; }
    
    std::string GetGridName() const { return _gridName; }
    std::vector<GridPoint<DataType>> GetPointArray() const { return _pointArray; }
    std::vector<CellType> GetCellArray() const { return _cellArray; }
};
